﻿using LxBasic;
using System.Diagnostics;
using System.Security.Cryptography;
using System.Xml;

namespace LxCore
{
    public partial class NSEquation : EquationWithDiffTerm
    {
        protected ConvectionSchemes ConvectionScheme;
        public IField<double> Conv = ConstField.Double1;//流动项前面的系数，动量方程为1，对于能量方程为比热容Cp     
        protected double[] MassImbalance;
        protected double[] InnMassFlux;
        protected double[] BouMassFlux;

        #region 设置矩阵
        public void SetZone_ConvTerm(IField<double> Conv, ConvectionSchemes convectionScheme, MassFluxBase massBase)
        {
            this.Conv = Conv;
            this.ConvectionScheme = convectionScheme;
            this.InnMassFlux = massBase.InnMassFlux;
            this.BouMassFlux = massBase.BouMassFlux;
            this.MassImbalance = massBase.MassImbalance;
        }
        public void SetBound_FlowValue(IField<double> baseValue, FvBound bound)
        {
            CheckCreator.Instance(ref EquationBounds[bound.ID], out FlowValueBound b);
            b.SetFlowValue(bound, baseValue);
        }
        #endregion
        #region 供外部调用的计算过程
        public override void FormABOnFaces(ThreadInfo threadInfo)
        {
            //计算内部面 A_LU
            FormAbOnFaces_inner(threadInfo);
            //计算边界面的
            for (int i1 = 0; i1 < EquationBounds.Length; i1++)
            {
                switch (EquationBounds[i1])
                {
                    case FluxWallBound fluxbou: formABOnFluxWallBound(Zone.Bounds[i1], fluxbou, threadInfo); break;
                    case FixValueWallBound fixv: formABOnFixValueWallBound(Zone.Bounds[i1], fixv, threadInfo); break;
                    case FlowValueBound flow: formAbOnFlowBound(Zone.Bounds[i1], flow, threadInfo); break;
                }
            }
        }

        public override void FormABOnCells(ThreadInfo threadInfo)
        {
            //massImbalance要加到系数上，例如流出的多，需要虚拟流入
            (int i1, int endp1) = Zone.Inner.CellIdx.Split(threadInfo);
            if (Conv.IsParameter(out double d) && d == 1)
            {
                if (Source_a.IsEqual0())
                {
                    if (Source_b.IsEqual0())
                    {
                        for (; i1 < endp1; i1++)
                        {
                            ref FvInnCell cell = ref InnCells[i1];
                            double massimb = MassImbalance[i1];
                            if (double.IsNegative(massimb))
                            {
                                ADiag[i1] = CalAdiag(ref cell) - massimb;
                                B[i1] = cell.CalIfOwnerAddElseMunus(InnFaceShare, BouFaceShare);
                            }
                            else
                            {
                                ADiag[i1] = CalAdiag(ref cell);
                                B[i1] = cell.CalIfOwnerAddElseMunus(InnFaceShare, BouFaceShare) + massimb * InnX[i1];
                            }
                        }
                    }
                    else
                    {
                        var sou_b = Source_b.GetCellVarArr(Zone.Inner);
                        for (; i1 < endp1; i1++)
                        {
                            ref FvInnCell cell = ref InnCells[i1];
                            double massimb = MassImbalance[i1];
                            if (double.IsNegative(massimb))
                            {
                                ADiag[i1] = CalAdiag(ref cell) - massimb;
                                B[i1] = cell.CalIfOwnerAddElseMunus(InnFaceShare, BouFaceShare) + sou_b.GetCellVar(i1) * cell.Volumn;
                            }
                            else
                            {
                                ADiag[i1] = CalAdiag(ref cell);
                                B[i1] = cell.CalIfOwnerAddElseMunus(InnFaceShare, BouFaceShare) + massimb * InnX[i1] + sou_b.GetCellVar(i1) * cell.Volumn;
                            }
                        }
                    }
                }
                else
                {
                    var sou_b = Source_b.GetCellVarArr(Zone.Inner);
                    var sou_a = Source_a.GetCellVarArr(Zone.Inner);
                    for (; i1 < endp1; i1++)
                    {
                        ref FvInnCell cell = ref InnCells[i1];
                        double massimb = MassImbalance[i1];
                        if (double.IsNegative(massimb))
                        {
                            ADiag[i1] = CalAdiag(ref cell) - massimb - sou_a.GetCellVar(i1) * cell.Volumn;
                            B[i1] = cell.CalIfOwnerAddElseMunus(InnFaceShare, BouFaceShare) + sou_b.GetCellVar(i1) * cell.Volumn;
                        }
                        else
                        {
                            ADiag[i1] = CalAdiag(ref cell) - sou_a.GetCellVar(i1) * cell.Volumn;
                            B[i1] = cell.CalIfOwnerAddElseMunus(InnFaceShare, BouFaceShare) + massimb * InnX[i1] + sou_b.GetCellVar(i1) * cell.Volumn;
                        }
                    }
                }
            }
            else
            {
                var conv = Conv.GetCellVarArr(Zone.Inner);
                if (Source_a.IsEqual0())
                {
                    if (Source_b.IsEqual0())
                    {
                        for (; i1 < endp1; i1++)
                        {
                            ref FvInnCell cell = ref InnCells[i1];
                            double massimb = MassImbalance[i1] * conv.GetCellVar(i1);
                            if (double.IsNegative(massimb))
                            {
                                ADiag[i1] = CalAdiag(ref cell) - massimb;
                                B[i1] = cell.CalIfOwnerAddElseMunus(InnFaceShare, BouFaceShare);
                            }
                            else
                            {
                                ADiag[i1] = CalAdiag(ref cell);
                                B[i1] = cell.CalIfOwnerAddElseMunus(InnFaceShare, BouFaceShare) + massimb * InnX[i1];
                            }
                        }
                    }
                    else
                    {
                        var sou_b = Source_b.GetCellVarArr(Zone.Inner);
                        for (; i1 < endp1; i1++)
                        {
                            ref FvInnCell cell = ref InnCells[i1];
                            double massimb = MassImbalance[i1] * conv.GetCellVar(i1);
                            if (double.IsNegative(massimb))
                            {
                                ADiag[i1] = CalAdiag(ref cell) - massimb;
                                B[i1] = cell.CalIfOwnerAddElseMunus(InnFaceShare, BouFaceShare) + sou_b.GetCellVar(i1) * cell.Volumn;
                            }
                            else
                            {
                                ADiag[i1] = CalAdiag(ref cell);
                                B[i1] = cell.CalIfOwnerAddElseMunus(InnFaceShare, BouFaceShare) + massimb * InnX[i1] + sou_b.GetCellVar(i1) * cell.Volumn;
                            }
                        }
                    }
                }
                else
                {
                    var sou_b = Source_b.GetCellVarArr(Zone.Inner);
                    var sou_a = Source_a.GetCellVarArr(Zone.Inner);
                    for (; i1 < endp1; i1++)
                    {
                        ref FvInnCell cell = ref InnCells[i1];
                        double massimb = MassImbalance[i1] * conv.GetCellVar(i1);
                        if (double.IsNegative(massimb))
                        {
                            ADiag[i1] = CalAdiag(ref cell) - massimb - sou_a.GetCellVar(i1) * cell.Volumn;
                            B[i1] = cell.CalIfOwnerAddElseMunus(InnFaceShare, BouFaceShare) + sou_b.GetCellVar(i1) * cell.Volumn;
                        }
                        else
                        {
                            ADiag[i1] = CalAdiag(ref cell) - sou_a.GetCellVar(i1) * cell.Volumn;
                            B[i1] = cell.CalIfOwnerAddElseMunus(InnFaceShare, BouFaceShare) + massimb * InnX[i1] + sou_b.GetCellVar(i1) * cell.Volumn;
                        }
                    }
                }
            }
            //ADiag[i1] = CalAdiag(ref cell) - massimb - sou_a.GetCellVar(i1) * cell.Volumn;
            //B[i1] = cell.CalIfOwnerAddElseMunus(InnFaceShare, BouFaceShare) + sou_b.GetCellVar(i1) * cell.Volumn;
        }
        public override void UpdateBoundsValue(ThreadInfo threadInfo)
        {
            for (int i1 = 0; i1 < EquationBounds.Length; i1++)
            {
                switch (EquationBounds[i1])
                {
                    case FixValueWallBound fixValue: SetXs(fixValue.FixValue, Zone.Bounds[i1], threadInfo); break;
                    case FluxWallBound fluxbou: updateFluxBound(Zone.Bounds[i1], fluxbou, threadInfo); break;
                    case FlowValueBound flow: updateFlowBoundValue(Zone.Bounds[i1], flow, threadInfo); break;
                }
            }
        }
        #endregion
        #region 内部被调用的计算方法
        protected void FormAbOnFaces_inner(ThreadInfo threadInfo)
        {
            var gammaf = Gamma.GetFaceVarArr(Zone.Inner);
            (int i1, int endp1) = Zone.Inner.FaceIdx.Split(threadInfo);
            if (Conv.IsParameter(out double d) && d == 1)
            {
                switch (ConvectionScheme)
                {
                    case ConvectionSchemes.FirstOrder:
                    default:
                        for (; i1 < endp1; i1++)
                        {
                            ref FvInnFace face = ref InnFaces[i1];
                            double massFlux = InnMassFlux[i1];
                            double gamma = gammaf.GetFaceVar(i1);
                            if (double.IsNegative(massFlux))
                            {
                                SetALU(i1, ref face, gamma, massFlux, 0, 1, 0);
                            }
                            else// massflux >= 0
                            {
                                SetALU(i1, ref face, gamma, massFlux, 1, 0, 0);
                            }
                        }
                        break;
                    case ConvectionSchemes.SecondOrder:
                        for (; i1 < endp1; i1++)
                        {
                            ref FvInnFace face = ref InnFaces[i1];
                            double massFlux = InnMassFlux[i1];
                            double gamma = gammaf.GetFaceVar(i1);
                            if (double.IsNegative(massFlux))
                            {
                                var f = GetFactor_FluxNegative(i1, ref face).CalFactor_SOU_CD_DC();
                                SetALU(i1, ref face, gamma, massFlux, f.fD, f.fC, f.fConst);
                            }
                            else// massflux >= 0
                            {
                                var f = GetFactor_FluxPositive(i1, ref face).CalFactor_SOU_CD_DC();
                                SetALU(i1, ref face, gamma, massFlux, f.fC, f.fD, f.fConst);
                            }
                        }
                        break;
                    case ConvectionSchemes.MUSCL:
                        for (; i1 < endp1; i1++)
                        {
                            ref FvInnFace face = ref InnFaces[i1];
                            double massFlux = InnMassFlux[i1];
                            double gamma = gammaf.GetFaceVar(i1);
                            if (double.IsNegative(massFlux))
                            {
                                var f = GetFactor_FluxNegative(i1, ref face).CalFactor_MUSCL_DC();
                                SetALU(i1, ref face, gamma, massFlux, f.fD, f.fC, f.fConst);
                            }
                            else// massflux >= 0
                            {
                                var f = GetFactor_FluxPositive(i1, ref face).CalFactor_MUSCL_DC();
                                SetALU(i1, ref face, gamma, massFlux, f.fC, f.fD, f.fConst);
                            }
                        }
                        break;
                }
            }
            else
            {
                var convf = Conv.GetFaceVarArr(Zone.Inner);
                switch (ConvectionScheme)
                {
                    case ConvectionSchemes.FirstOrder:
                    default:
                        for (; i1 < endp1; i1++)
                        {
                            ref FvInnFace face = ref InnFaces[i1];
                            double massFlux = InnMassFlux[i1] * convf.GetFaceVar(i1);
                            double gamma = gammaf.GetFaceVar(i1);
                            if (double.IsNegative(massFlux))
                            {
                                SetALU(i1, ref face, gamma, massFlux, 0, 1, 0);
                            }
                            else// massflux >= 0
                            {
                                SetALU(i1, ref face, gamma, massFlux, 1, 0, 0);
                            }
                        }
                        break;
                    case ConvectionSchemes.SecondOrder:
                        for (; i1 < endp1; i1++)
                        {
                            ref FvInnFace face = ref InnFaces[i1];
                            double massFlux = InnMassFlux[i1] * convf.GetFaceVar(i1);
                            double gamma = gammaf.GetFaceVar(i1);
                            if (double.IsNegative(massFlux))
                            {
                                var f = GetFactor_FluxNegative(i1, ref face).CalFactor_SOU_CD_DC();
                                SetALU(i1, ref face, gamma, massFlux, f.fD, f.fC, f.fConst);
                            }
                            else// massflux >= 0
                            {
                                var f = GetFactor_FluxPositive(i1, ref face).CalFactor_SOU_CD_DC();
                                SetALU(i1, ref face, gamma, massFlux, f.fC, f.fD, f.fConst);
                            }
                        }
                        break;
                    case ConvectionSchemes.MUSCL:
                        for (; i1 < endp1; i1++)
                        {
                            ref FvInnFace face = ref InnFaces[i1];
                            double massFlux = InnMassFlux[i1] * convf.GetFaceVar(i1);
                            double gamma = gammaf.GetFaceVar(i1);
                            if (double.IsNegative(massFlux))
                            {
                                var f = GetFactor_FluxNegative(i1, ref face).CalFactor_MUSCL_DC();
                                SetALU(i1, ref face, gamma, massFlux, f.fD, f.fC, f.fConst);
                            }
                            else// massflux >= 0
                            {
                                var f = GetFactor_FluxPositive(i1, ref face).CalFactor_MUSCL_DC();
                                SetALU(i1, ref face, gamma, massFlux, f.fC, f.fD, f.fConst);
                            }
                        }
                        break;
                }
            }
        }
        protected ConvectionDiscretization GetFactor_FluxPositive(int i1, ref FvInnFace face)
        {
            return new ConvectionDiscretization(InnX[face.OSideCell], InnX[face.NSideCell], 2 * Vector3d.DotProduct(Grad[face.OSideCell], face.OToN), face.gOwner);
        }
        protected ConvectionDiscretization GetFactor_FluxNegative(int i1, ref FvInnFace face)
        {
            return new ConvectionDiscretization(InnX[face.NSideCell], InnX[face.OSideCell], (-2) * Vector3d.DotProduct(Grad[face.NSideCell], face.OToN), 1 - face.gOwner);
        }
        protected void SetALU(int i1, ref FvInnFace face, double gamma, double massFlux, double a, double b, double c)
        {
            double _gamadiff = -gamma * face.gDiff;
            AInnO[i1] = _gamadiff + massFlux * b;
            AInnN[i1] = _gamadiff - massFlux * a;
            InnFaceShare[i1] = gamma * face.CalNonOrthogonalityCorrectionOfInner(Grad, InnX) - massFlux * c;
        }
        protected void formAbOnFlowBound(FvBound bound, FlowValueBound flowBound, ThreadInfo threadInfo)
        {
            (int i1, int endp1) = bound.FaceIdx.Split(threadInfo);
            var conv = Conv.GetFaceVarArr(bound);
            var gf = Gamma.GetFaceVarArr(bound);
            var factors = flowBound.BouFactors;
            int start = bound.FaceIdx.Start;
            for (; i1 < endp1; i1++)
            {
                ref var face = ref BouFaces[i1];
                double massflux = BouMassFlux[i1] * conv.GetFaceVar(i1);
                if (massflux > 0)//由明确的流出时，确定为流出边界
                {
                    //向外流出，边界条件为：梯度为0，因此只考虑流动，不考虑扩散
                    double innX_BouX = InnX[face.OSideCell] - BouX[i1];
                    ABouH[i1] = -massflux;
                    BouFaceShare[i1] = massflux * innX_BouX;// (InnX[face.OSideCell] - BouX[i1]);//用downwind格式
                    factors[i1 - start] = -(Vector3d.DotProduct(Grad[face.OSideCell], face.FaceNOrth) - innX_BouX * (face.gT - face.gDiff))
                        / face.gDiff;
                }
                else
                {
                    //流入，边界条件为：固定值，同时考虑流动和扩散
                    double gamma = gf.GetFaceVar(i1);
                    double gammadiff = gamma * face.gDiff;
                    double xbou = BouX[i1];
                    double nonOrthCor = face.CalNonOrthogonalityCorrectionOfBound(Grad, InnX, xbou);
                    BouFaceShare[i1] = gamma * nonOrthCor - (massflux - gammadiff) * xbou;
                    ABouH[i1] = -gammadiff;
                    factors[i1 - start] = -nonOrthCor / face.gDiff;
                }
            }
        }
        protected void updateFlowBoundValue(FvBound bound, FlowValueBound flowBound, ThreadInfo threadInfo)
        {
            (int i1, int endp1) = bound.FaceIdx.Split(threadInfo);
            var bf = flowBound.FlowValue.GetCellVarArr(bound);
            var factors = flowBound.BouFactors;
            int start = bound.FaceIdx.Start;
            for (; i1 < endp1; i1++)
            {
                if (BouMassFlux[i1] > 0)//由明确的流出时，确定为流出边界
                {
                    BouX[i1] = InnX[BouFaces[i1].OSideCell] + factors[i1 - start];
                }
                else
                {
                    BouX[i1] = bf.GetCellVar(i1);
                }
            }

        }
        #endregion
        public class FlowValueBound : EquationBoundBase
        {
            public IField<double> FlowValue;
            public double[] BouFactors;
            public void SetFlowValue(FvBound bound, IField<double> fixValue)
            {
                CheckCreator.ArrayStrict(ref BouFactors, bound.FaceIdx.GetCount());
                this.FlowValue = fixValue;
            }
            public override void Initialize(EquationBase inEquation, FvBound inBound, ThreadInfo threadInfo)
            {
                inEquation.SetXs(FlowValue, inBound, threadInfo);
            }
        }

    }
}
